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ABSTRACT 


The most effective way to design VLSI MOSFET structures is to 
use a sophisticated complex two dimensional model. Such a model, 
using the basic semiconductor equations along with the accurate 
models for the physical parameters of the basic equations, has 
been described here. Numerical simulation of this MOSFET model 
using finite difference method has also been described. Low 
computation cost algorithms have been used to develop an 
easy-to-use software package for planar MOSFET simulation. 
Simulated results are presented for few types of short channel 
MOSFETs to show the power of the device simulation to predict the 
behaviour of the device. Finally, the modifications of this MOSFET 
simulation model and the possibility to integrate this device 
simulation with the process and the circuit simulations have been 


proposed . 
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INTRODUCTION 

The development of new semiconductor technologies and novel 
semiconductor device structures have traditionally been guided by 
an experimental approach . Starting from an established process 
sequence , fabrication steps are changed together with geometrical 
dimensions ( feature size ). The modified process is then realized 
by fabricating several lots . The finished devices are tested to 
insure their performance conforms to the initial design . This 
approach usually includes several iterations of the process to 
testing loop . Uith the advent of increasingly complex integrated 
circuits , the traditional empirical approach becomes expensive 
and time consuming . 

There are analytical models for the device and the process 
steps available as an alternative approach towards the design and 
development of integrated circuits . But all of them use certain 
assumptions and regional approximations , which are often so 
restrictive that the predictability of the performance of the 
devices in the VLSI complexity has turned out to be insufficient. 

In order to characterize the modern devices in a reasonable 
way , higher order n.xmMX'ica.l nuadtals emerge out to be the only 
solution . The extraction of the solution of these numerical 
models , which approximate the actual physics going on in the 
devices , in a fabrication process step or in the fabricated 
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circuit closely, is popularly called as . When 

simulation is done for a device , it is called stm-ulation ; 

similarly for the fabrication process and for the circuit 
performance are' called px-oc&ss arui cix-cuit stmiulaticms 
respectively . 

The device simulation when coupled with the process 
simulation and circuit simulation , makes design and development 
process of VLSI circuits very cost and time effective compared to 
the experimental approach described before . Basically , the 
simulation softwares are guided by the principles of robustness , 
speed of solution , low computation cost and easy user access 

This thesis work deals with the development of a device 
simulation software and the planar Metal Oxidgt Pi«ldL E//»ct 
Tr-cinst&tor- ( MOSFET ) has been chosen as the device for 
simulation, because it is the semiconductor device widely used in 
the VLSI circuits . 

The dJewiee mi.mula.tion. basically means the solution of three 
coupled non-linear partial differential equations, i.e. Poisson’s 
equation and two current continuity equations , inside the domain 
of' the device ( here MOSFET ) with proper boundary conditions and 
with proper models of the physical parameters appeared in the 
equations mentioned . Basic partial differential equations (called 
the basic semiconductor equations) which describe the behaviour of 
the device and their boundary conditions are dealt in ofuxptmx^ 2 
This chapter also deals with the scaling of these equations for 
better numerical stability and accuracy . 
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C?i£tpt&r- 3 contains the models for the physical parameters 
like carrier mobilities, doping profile etc, and the assumptions 
imposed on the numerical model of MOSFET to make the computer 
solution faster without major loss of accuracy . 

Numerical methodologies to solve the basic coupled 
semiconductor equations ( i.e., d.iscr■^stista.tion and m&sh 
g&nsra.tion.') are given in chapter 4 . The algorithm for fast and 
low cost solution methods and post-processing of the results are 
dealt in chapterr- 5. 

Few example devices simulated by the device simulation 
software developed in this thesis work are described in the 6^ 
chapter . 

The concluding chapter ( chapter 7) outlines the scope of the 
extension and the modification of the simulation software . An 
apperwifx is also included to explain how to use this device 
simulation software . 
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MATHEMATICAL FORMULATION 

2. 1 . Basic Equations : - 

In order to accurately analyze any semiconductor structure, 
taking all the effects (non-ideal) into account, the basic 
semiconductor equations must be solved. The following are the 

equations governing inside semiconductor : 

2.1a. Poisson's aquation i 

2 ^ 

^ W = ~ - C P - n + N^-N^) (2.1) 

where, yf : electrostatic potential 

p,n : hole, electron concentrations 
N. -N ; total ionized Impurity 

0. fit 

« : permittivity of semiconductor 


2.1b. Continuity aquations 1 


1 

q 


V. J 


n ” Wt 


+ R 


1 

q 


v. j 

p 



( 2 . 2 ) 


R 


(2.3) 
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where, J ,J : electron,hole current density 
’ n’ p ’ ^ 

R : recombination rate 

q : electronic charge 

2.1c. Current, relations : 

= -q ( n ) (2.4) 

= -q C /Up p VV i ) (2.5) 

where, /u ,p : electron, hole mobilities 
’ '^n '^p ’ 

D_^ : electron, hole diffusion coefficients 
n, p 

:both of them are functions of electric field, dopant & carrier 
concentrations . 

Under non-degenerate condition D & P are related by 

p , n P > n 

Einstein relation : 

( 2 . 6 ) 

where, k :Boltzmann's constant & T :absolute temperature 
Under non-degenerate condition n & p are related to v' by : 

q(W-<^_)/fcT q(,^ -V)/kT 

n= e & p= nj^ e ^ (2.7) 

where, <6^, 4 ^^ : quasi fermi potentials for electron & hole, 

2 

under equilibrium condition & np= n^ . Uhere n^^ is 

intrinsic carrier concentration. 
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Nov returning to (2.4) & (2.5) first parts of r.h.s. are 
drift and second parts are diffusion components of currents. 
Combining these two parts (i.e., operating '7 over n & p and using 
(2.6) & (2.7) ) we get. 


V/V_ -^/^T 

J = q {D n, e ^ V (e ^) } 

n ^ ' n i ^ * 


-q{Dp n. e " V (e'P ^) ) 




( 2 . 8 ) 


(2.9) 


where, VY=kT/q : thermal voltage 

Now combining (2. 8), (2. 9), (2.2) & (2.3) to get combined 

equation for steady state [i.e., d/dts 0 ] 


V. {D^ V(e 


) ) = R 


( 2 . 10 ) 


-W/V 4, /V 

V.{Dp n^e " V(e P ^) ) = r (2.11) 

These (2.1), (2.10), (2.11) are the basic equations for the 
numerical model of semiconductor devices. 


2.2. MC»»MALIZATIC»t i> 

In order to get accurate computer solution, we need to scale 
the basic equations in such a way, that the dependent variables 
and independent variables vary within allowable limits and become 
dimensionless. One such normalisation scheme is proposed by De 
Marl [2] is adopted and the scaling factors are tabulated in the 


tabl6-2 . 1 . 




Scalinfi factors 
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Normalizing with the scaling factors (2.1), (2.10) & (2.11) 

become : 

= (n - p - C) (2.12) 

+ 

where,— C=N -N. , n= e & 

' ad' 

W 

V.{r^ e V(6 '^ )} = R 

4r 

'^■irp e v(6 ‘^ )) = R 

where, :normalized mobilities of electron & hole 

All th& variaSbl&s and o/ th* main. &q.\iatiorts wx^itt&n 

from (2.12) caud omaax-ds otre rvox-mattse-dL. Separate no tat tone ore not 
tijsec? /or t?ie normalfsed «giuantftfee to distinguish, thsm. /rom. the 
itneoaled equations t for- conusnisncs of -iax-iting. 


p- e 


C^p-V) 


(2.13) 


(2.14) 


Z. 3. CHOICE OF DEPENIHEINT VARIABLES t - 


For analytical purpose it is often helpful to use dependent 
variables other than n, p in the basic equations (2.12) through 
(2.14). One set, very frequently used, is {V.u.v} 
[SlotboomCB], Adachi et altll] ,Hiemier[12] & others]; which relates 
to the original set (Yfti.p) by : 

n*e .u & u=e 
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-w 4 > 

P=6 .V & V=6 

The major advantage is, it makes the continuity equations 
linear self-adjoint elliptic partial differential equations 
(P.D.E.)- Rugged solution methods are available for this type of 
problem. It is not used here as the major drawback of these 
variables is that, u and v need enormous dynamic range of 
numbers to represent them on computer [u (v) varies more than 33 
orders of magnitude for variation in -1 to +1 range ].So 
this set is limited to low voltage applications to semiconductor 
devices . 

Another set of dependent variables is , used as 
natural consequence from (2.12) through (2.14). An advantage here 
is that all the dependent variables are of same order of 
magnitude. Still there is a major disadvantage, all the equations 
become non-linear. This set has been used by Gummel[l], 
Heinerzhagen et al [16] and others [19]. 

The most suitable set is (v',n,p}. The advantages are same as 
the set {v'fWfV}; only poisson’s equation is rtton-linoar elliptic 
P.D.E. , and continuity equations are serif asLj'oirtt olliptic 
P.D.E. . Moreover the dynamic range for n and p is smaller [varies 
16 orders of magnitude at the maximum ] than the set of u and v. 
Fitchner ,Rose & Bank [8] pointed out the chance of getting 
inadmissible and unphysical negative carrier concentrations as the 
solution of basic equations due to discretization error in 
difference scheme [discussed later]. But the proper grid structure 
and proper discretization scheme for continuity equation eliminate 
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this discrepancy totally. This set was used by Van Dorp6[4], De 
Mari [ 2 ] ,Selberherr et al[13]. have chosen this set of 

dependent variables {y/^n^p}. 

Mock[6], Toyabe[143, Yamaguchi [ 20 ] and others used stream 
functions as dependent variables in continuity equations. But in 
this formulation there is no provision for inclusion of 
recombination-generation term. 


2.4. BOUNDARY CONDITIONS i- 

Before imposing boundary conditions, we need to discuss the 
structure of the device and the coordinate system used. Referring 
to figure(2.1) ’x’ is away from Si-SiO^ interface, and ’y’ is from 
source to drain contact. Since there is no change in n, p & yf in 
the bulk of semiconductor when MOSFET operates in subthreshold, 
linear and active mode computer simulation is done upto a 
reasonable depth ( 10 to 20 times the depletion width ). 

The boundaries present in the MOSFET structure [fig. (2.1)] 
can be divided into two parts : 

2. 4a. Raral boundarlas i 

In fig(2.1) ab, Ij, be, cd & kl fall under this category. Now 
the boundary condition for dependent variables are set as follows: 
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(TilJ Ptste-ntial .--At ohmic contacts [ab rsotnr-ce , ij :gate, 
cd:dr-ain, kl:i>£>dty] we know the applied potentials. So the value of 
yf at these places (Dirichlet boundary ), are given by. 


¥ 


V 


applied 


+ V 


bi 


(2.15) 


Where, : built-in junction potential ( applied bulk 
voltage is taken as reference, so is added with source and 
drain voltages). 

At Si-Si02 interface (be) Gauss’s law is applicable. 


^ c- fE ^ = 


^ 55E ^ 


+ G 


tni 


(2.16) 


where, Q. : interface trapped charge assumed to be present 

LITit 

at the interface as delta function. 

4£> . & £ : are permittivities of Silicon and 

at ox 

Si licon-di -Oxide 
- ^ : is electric field 


CiiJ} El&ctr-on & Wola - As proposed by H. K. Gummel[13, we 
took neutral space charge density and thermal equilibrium for 
carriers at ohmic contacts, i.e.. 


C 



n - p -C = 0 ; 


& n.p = 1 
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Solving these equations : 



n = / + 1 j + CC/2) 

(2.17) 

- (C/2) 

It is also clear from np = 1 that, are same at the 

ohmic contacts. 

At the interface, there is no component of current density 
normal to the boundary (he) i.e., 

j =0 «« ^px=° 

nx 

Corresponding boundary conditions for n & p are derived 
f om (2 18) during discretizat ion of continuity equations. 

2. 4,b. Fictitious boundaries i 

They are not the physical limits of the device, but imposed 
as the boundary of simulation area [aek, dfl, ib, ic, ef of 

fig. (2-1) ) - 

At these boundaries no changed in dependent variables 

flowing oat of 

these boundaries. 
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0 , 



Sti 
dy " 


0 


At the vertical fictitious boundaries. 


( 2 . 19 ) 




0 , 




( 2 . 20 ) 


At the horizontal fictitious boundaries. 

Now, with these scaled basic equations and boundary 
conditions in one hand, we are going to establish the practical 
numerical model for the MOSFET and other physical parameters like, 
mobility, doping profile. Intrinsic carrier concentration etc. in 


the next chapter. 
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MODELS & ASSUMPTIONS 

This chapter consists of the suitable numerical models 
for the MOSFET and the physical parameters, i.e., carrier 
mobilities, intrinsic carrier concentration, doping profile of 
the device under consideration etc. 

3. 1 . NUMERICAL MODEL FOR MOSFET t - 


The model hierarchy consists of three different steady 
state MOS models . The highest level ( level 3 ) of hierarchy 
requires numerical solution of two dimensional Poisson’s equation 
and both of the continuity equations in two dimension with 
appropriate boundary conditions and physical models . Only level 3 
model can handle avalanche generation and breakdown in MOSFET. 
Uhen avalanche generation starts, the coupling between Poisson’s 
and continuity equations becomes very strong; so all the basic 
equations are to be solved simultaneously ( Meinershagen et al 
[16] ). This simultaneous solution is very expensive in terms of 
computer memory and C.P.U time. The two lower level models are 
valid when there is no generation-recombination. 

Level 2 model assumes that there is no minority carrier 
current [This assumption is supported Selberherr et al :MINIM0S1 
[13], Toyabe et al , CADDET [17], De la Moneda [19] and others ]. 
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For n-channel MOSFET assumption is : 

J = 0 in the entire device. 

P 

So, level 2 model simplifies down to solution of 2-d 
Poisson's and 2-d majority carrier (electron for n-channel) 
continuity equation . 

Level 1 is the simplest model, assumes no minority 
carrier current and majority carrier current flows only in lateral 
(along source to drain ) direction . So : 

n 

^ =0 in the channel region. 

With this simplification, solution work tackles only 2-d 
Poisson’s equation and 1-d majority carrier continuity equation. 
This is quite accurate model in subthreshold and linear regions of 
MOS I-V characteristics. 

We have used, le-veL 2 mod&l , as our main interest is to 
simulate MOS in ^xtblhr^s^holat, lin&ctr- and activa zone . Moreover 
level 2 model is tfuita- ^aonomic in terms of camput^x- jnemory and 
C.P.U. timie . This model is claimed to be reasonably close with 
the experimental results by several authors [13], [14], [19], 

[ 20 ]. 

The level 2 model used in this simulation is written 
again: ^ 

<r<aL> N— channel MOS : - 


s=(n-p“C) 


(3.1a) 
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V.3 = 0 
n 

(3.1b) 

J = 0 

(3.1c) 

P 

i.e., ^ = constant 

^P 

(3. Id) 

Ct>y P— channel MOS - 


2 

7V' = (n-p-C ) 

( 3 . 2a) 

7.J = 0 

P 

(3.2b) 

Jn = 0 

(3.2c) 

i.e., df = constant 

(3. 2d) 


3.2. OTHER AlSEUMPTIONS GENERALLY USrED WITH LEVEL 2 MODEL s- 

ay of p&r-mi 1 1 ivi ty : 

e .= constant ,« , = constant C3>3a) 

sv ®*-°2 

cay Total ion.i£iation. o/ otl I tho imfuMX'i t i&s 

C = = n/ - n/‘ (3.3b) 

day No han/I gOLp nax'x^owing : 

19 

Since doping crosses 10 /cc only at positions 
very close to the ion implanted source and drain contacts, other 

19 

parts have much less than 10 /cc impurity concentration. So we 
assume, 

AE =*0 i.e.,An, = 

& i 


0 


(3.3c) 
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Homo^&njet ty of t&mp&x-ot'ux'O : 

T = constant over the whole body of the device 

(3. 3d) 


Cvy No T^-ocomhinjatiorL ctrud. gonei-otion. : 

This forms the validity of the level 2 model for MOSFET. 
R = 0 C3.3e) 

3.3 MODELLING <W THE PHYSICAL PARAMETERS i- 

3. 3a. Thermal voltage i 



-23 * 

where, k : Boltzmann’s constant =1.381 xlO Joule/ K 
T : Absolute temperature ( •■K) 

-19 

q : electronic charge = 1.6 • 10 coulomb . 


3. 3b. Intrinsic carrier concentration C n^ 3 i 

n^CT) = 3.88 k 10^^ T 6 ( / c c ) 

3.3c. Carrier Mobility i 


(3.5) 


Accurate mobility values are required for the purpose of 
predictive simulation because of the multiplicative dependence of 
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current upon mobility . Various mechanisms influence the carrier 
mobility ; they are scattering due to thermal vibration of 
lattice, ionized impurities, scattering of electron and hole, 
velocity saturation due to parallel electric field and the effect 
of transverse electric field etc . 

There are several mobility models derived from 
theoretical background, approximately taking all the effects and 
their interactions . From simulation point of view, 
models fit to the experimental data are of more importance . These 
type of models are used by several authors [ Toyabe[17], 
Yamaguchi [20],Navon [21] ] . 

For the thisr-mal vibr-ation. of lattico ^ mobility 
degradation is taken as a simple power law [ tested and fit by 
Arora et al [22], Grove [23] ] : 

C T / 300“K ) 

Jr 

o 2 

where, = 1350 cm /v/s = 2,5 

u 480 cm^/v/s « = 2.5 

P P 

* The values are taken from Grove [23], similar expressions are 
supported in MINIMOS 1 [13] . 

The empirical model including the effect of 


(3.6a) 

(3.6b) 


<r£J» imptix-ity s>ca.t tor-ing 
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velocity lion. ot£ hi^h. olootr-ilo 

Cpaz-a.llerl to otiT-r-ont c^ia^r•ont flaw & 

<r £ £ £ J> sMX-faaa scat t&r-ing 

was proposed by Yamaguchi [20], stated as : 
A^CNg,Ep,E^) — fj . f(Ng,Ep) .gCE^) (3.7) 


f ialdL 


where, : impurity concentration, 

Ep : parallel electric field 

E^ : transverse electric field, often referred as 
’gate field’ . 

The term ’f’ represent the bulk mobility reduction 
factor . This term was proposed by Scharfetter and Gummel [3] : 


fCNg.Ep) =. [ 


N. 


1 + 


B 


(Ep/A) 


(Ng/S) + N 


(Ep/A) + F 


+ (Ep/B)^ j 


- 1/2 


(3.8) 


The first two terms account for impurity scattering 
mobility reduction. The constants S & N are explained in 
figureCS.l) . ’N’ is the 'corner' impurity concentration where 
mobility reduction begins. VS is the maximum reduction of log(p). 
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The values given by Schar fetter and Gummel [3] are : 



S 

N (/cc) 

Electron 

350 

3.10^^ 

Hole 

81 

4.ld*^ 

Li — 


* Table 3.1 * 


The explanation of the last terms are given [Thornber 
[24]] by rearranging (3.8), 


fCNfi.Ep) 


[ 


N 


1 + 


B 


(Ng/S) + N 


+ (iu^) 


(Ep/Cju^.A))^ 


{jt . E 


/(p'^.A) + F 


+ 



2 1 


(3.9) 


The term (/j^.B) represents the saturation velocity 
which is independent of impurity Impurity scattering . The term 
(p'^.A) can be identified as acoustic phonon velocity, responsible 
for warm carriers. Here ’F’ does not have any physical meaning, 
but a fitting factor of the empirical formula to the experimental 
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values of the constants at T=300‘K [33is 

tabl6(3. 2) . 



ACV/cw j 

' 

F 

B(v/cw) 

Electron 

3 . 5x 1 0 ^ 

8.8 

7.4 xlO^ 

Hole 

6 . Ik 1 0 ^ 

1 . 6 

2.5h10^ 


data . 

The 

tabulated in 


* Table 3.2 * 

The temperature dependence of B and A constants given by 
Selberherr [25], used here : 

7 -0.S7 

B = 10 (T/300*K) n electron (3. 10a) 

, -0 . 52 ^ 

B = 8.37k10® CT/300‘K) / hole (3.10b) 

Jr 

Temperature dependence of A is claimed [25] to be very 
weak, so its value in the table3.2 has been considered. 

The 'g' term is introduced in (3.7), so as to 
represent the mobility . The carrier motion is strongly 

affected by the roughness of the surface and the interface states. 
The surface mobility decreases with the increase in gate field 
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[Fang & Fowl6r[263 ]. An approximate formula which represents this 
property was proposed and established by Yamaguchi [20] . 

- 1/2 

gCE^) = C 1 + « ) C3.ll) 


The parameter a. is reciprocal of critical electric 
field .The value given by the same author as : 


<x 


1.539h 10 cm/v for electron 


5.35k10 


"5 CJ((s jV 


for hole 


In the model (3.7) carrier-carrier scattering has not 
been taken into account . In the program, this effect has been 
included by a very simple approach [ proposed by Engle and Dirks 
[27] ], 1.6., replacing by : 

Ng = 0.34 times doping concentration +0.66 times (n+p) 

(3.12) 

TTue o/ modisls huere 

&staiiblimh&€i hy ths axithor-s lih»B Ycmua^uchi t£01 » Sslb^r-hor-r- » 

Schao'f&tttfT- larud G^mawsl [3J and in atfmr- pap»x-tif. 

3.3d. Doping profllo i 

Accurate doping can only be calculated by process 
simulation programs like SUPREM [28], SUPRA [29] etc.; these are 
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the problems of process modelling and simulation . In this device 
simulation program, implanted Gaussian doping profile in the 
source, drain and channel implantations has been assumed . 

The vertical implanted Gaussian profile for source, 
drain and channel reads : 


2 

Dese (x-R ) 

C(x) = exp - .2 I (3.13) 

•*G AR P 

P 


2 

where. Dose : implantation dose in atoms/cm 

R : projected range of implanted profile in cm 

Jr 

ARp : straggle of the same in cm 


For source and drain lateral diffusion exact profile can only 
be given by process simulation . An approximate Gaussian profile 
has been assumed in lateral direction : 


Dose 


CCx,y)= 


exp 


, AR 
2n p 


(y-Lg) + (*-Rp) 


2 AR 


Dose 


exp 


<y-(Ly-L^)) + (x-R^) 


2 » 


2 AR 


^ ] 


(3.14) 
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where, L & L. : Source & drain contact lengths 


L : Lateral dimension of simulation area. 

y 

The values of R , AR are taken from graphs of 

P P 

quantities against implantation energies (in KeV) [25]. 


Recombination generation has not been modelled as it has 


these 


not 


been included in this simulation. 
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NUMERICAL TREATMENT 

The aysteia o£ baaic aetnlconductor equatlona cannot be aolved 
explicitly taking all the non-ideal effecta into account. 
Therefore, the solution must be calculated by means of numerical 
approaches. This chapter deals with these things. 

4,-1 BASIC CONCEPTS s- 


The numerical approach for the solution of differential or 
partial differential equations, consists of essentially three 
tasks. First one is, the domain, i.e., simulation geometry of the 
device, has to be partitioned into a finite number of subdomains, 
in which the solution can be approximated with a desired accuracy. 
Secondly, the differential equations have to beapproximat ed in 
each of the subdomains by algebraic equations which involve only 
values of the continuous dependent variables at discrete points. 
In this way, one obtains a fairly large system of algebraic 
equations with unknowns comprised of approximations of continuous 
dependent variables. The third task is to solve this system of 
equations . 

It is to be mentioned here that the exact analytical solution 
cannot be found by the method outlined above. Rather, we get an 
exact solution of the system of algebraic equations. If the 



subregions are small enough to ensure better fit to the exact 
values the numerical solution represents a good approximation to 
the solution of analytically formulated problems. 

There are many classical methods which propose the 
possibility of the partitioning of the domain (discretization). 
The are : 


(i) Finite Difference Method (F.D.M.) 

(ii) Finite Box Method (generalized F.D.M.) 

(iii) Finite element Method (F.E.M.). 

Ue used F.D.M. discretization as it is easier to tackle 
regular shaped geometry like that of MOSFET, than other methods. 


4.2 GRID FORMATION 

In F.D.M. the domain is discretized by straight lines (grids 
or mesh) parallel to X and Y axes. The spacing is kept 
non-uniform, closely spaced grids are used where there is a chance 
of the variation of any dependent variable to be stiff, and where 
the variation is less, grids are spaced widely apart. (More about 
the generation of non-uniform grid or mesh is discusses later.). 
One possible grid-structure is shown if Figure 4.1; it is used as 
the starting mesh in this program. There are NX lines in 
X-direction, NY in Y-direction, so there are (NX. NY) grid points, 
where we have to approximate the analytical differential 


equations . 






4.3 DISCRETIZATION i- 


Here, the differential equations are replaced by difference 
equations on the grid points by suitable discretization techniques 
[30], [31]. 

Ue used five point molecule (Figure 4.2) for the 
discretization of the 2-d partial differential equations (the 
basic equations discussed in the previous chapter). 

The notations used in the discretization are given as follows: 
Grid spacings : 

hi = ; i = 1 to NX - 1 (4.1a) 
ittj = y^^.^ - Yj ; j = 1 to NY - 1 (4.1b) 
Discrete point dependent variables : 


f^ j = f(x^, y^; i = 1 to NX, j = 1 to NY (4.2a) 


i + 172 , j 


*i**l+l 
= f ( 


, Yj); i * 1 to NX, j=l to NY (4.2b) 


^i, j+1/2 


^ C ^ > 


2 


); i 


1 to NX, j=l to NY (4.2c) 


Where, ’f' is any one of the dependent variables v, n or p; 
physical parameters C (doping density), or y^ (normalized 



mobility) or parameters like E (electric field), current density 


J or J . 
n p 


4. 3a Discretization of Poisson* s Equation s- 


In 2-D simulation, we assumed MOS with large width (W) and no 
variation of any of the dependent variables in Z-direction (normal 
to the paper). Hence the normalized Poisson’s equation becomes 
as : 


ay? 


a\ 


ay 


(e - e ** - C) 


(n - p - C) 


(4.3) 


It is linearized by Newton’s method [1], assuming. 


^ exact 


V + ot 


(4.4a) 


±ot 

0 


^ 1 ± ot 


(4.4b) 


Where, « is very small quantity, defined as the error in 
potential when we approximate (4.3) as linear one. 

With the assumptions (4.4a) and (4.4b), Poisson’s equation is 
linearized down to : 

£4 . = („ - p -C) - [£!| . 

aycr ay^ ax ay 


(4.5) 


The Poisson’s equation is solved iteratively by making a 
guess at K-th iteration for y/ , a is calculated from (4.5), the 
guess is refined as (4.6) and proceed to K+lst iteration. 


K+1 


V' 



+ a 


(4.6) 


K 

This is iterated until a becomes negligible. 

Now referring to Figure 4.2 of non-uniform 5 point molecule, 
y is expanded in Taylor’s series about the points (i+1/2, j) and 
(i-1/2, j) and subtracting these expansions, we get after 
rearrangement , 




h 


h 


where, 0(h) : of order h 


i-1 


O(h^) 


of order h 


Chf - 


hi_i)/24 


lET 


~hJZ^17T 


etc . 


(4.8) 


( 


Similarly, analogous expressions can be derived for 

55E'’i.j’ 

Using (4.7) and its analogues in (4.5) we get. 


W^i 



3 


r^i 

UxJ i + l/2, j“[dxj i-l/2, j . 

iiL..Wr^irTi-|.iiii '--"-' r- oim.-niiumr 


h , + h . , 
i 1-1 


0(h)(^) 

ax^ 


fa 

otl p 

a? 

+k 

yJi,j + l/2 U 

yji, j-l/2 




(n. .+p, .+C. .) 

X f i 1 , j 1 » j 


i + l/2, i-1/2 , j ^ 


^ ^i-1 


Sx 


3'i,} 


Xn m * Xtt * i4 A ^ f J 

j j-i 


(4.9) 


Sw 

Next. we derive aa equation (4.7) for ^^^i-fl/2 j 
etc , 


C— ) 

'‘dx''i + l/2, j 


C^) 

'•dx'’i + l/2, j 


+ o(h 


(4.10) 


and ao on. 

We aaauflie, h. and m . are amall enough to make the effect of 

* J 

0(h) or 0(m) times the third derivative of potential and similar 
higher order terms negligible [they are called the local 
truncation errors]. Combining (4.9) with (4.10) and its 

analogues, we get the difference equation for linearized and 
normalized Poisson’s equation [excluding local truncation errors] 


as follows : 



3 


a,., , - a, , ot, . - a, ^ . a. ... - a. . 


h 


i-1 


tn 


m . , 
J-1 


h. + h, , 
1 i-1 


®j “j-l 


■ "i.j ■*■ ^L,}^ = ~ Pi,j " ^i,j^ 


T i + l , j . j '‘'i.j ^i-l,j ^i,j + l ~^i,j ^i.j ~^i,j-l. 


tl, 

1 


^i-1 


m 


“j-i 


h. + h. , 
1 1-1 


m . + m . . 
J J-1 


C4.ll) 


Local truncation errors (L.T.E.) in the difference equation 
(4.11) are as : 

Truncation error in a (dropping subscripts) 


T ^ 0(h) 
a 


*3 

a 


Sx' 


+ 0(in) 


3 

S a 


Sy' 


(4.12) 


L.T.E. in ^ is just similar to (4.12), with « replaced by y. 

From equation (4.12) it is clear that, L.T.E. is linearly 
proportional to the grid spacing. So, closely spaced grid is to 
be set, where v' (so ") varies sharply to get better approximation 
to exact solution. Widely spaced grids are used where v varies 
very little (in the bulk of MOS) so as to decrease the number of 
grid points (to minimize C.P.U. time on computation). 

One important point is to be noted here; by comparing 
equation (4.5) and difference equation (4.11) one can conclude 



that v* (and « also) is taken to be varying linearly between two 
neighbouring grid points, during the derivation o£ difference 
equation. 

Now, we derive the difference equation for Poisson’s equation 
at the boundaries. 


(i) At the ohroic contacts of source, drain and gate, v* is 
known as 


V = V + V 

applied built 

So is known and from 

oi = 0 

Setting this condition 
equation at ohmic contacts. 


in 


(4.12) 


equation (4.4a), it is clear that 

(4.13) 

in equation (4.11) we get difference 


(ii) At the fictitious boundaries [j = 1 & NY and = 
simulation depth], we assumed in chapter 2 that no electric field 
is going out of the device : 


^ da 
dy dy 


0 at j 


1 & NY for 1 s: i NX 


(4.14a) 


Replacing it by difference approximation. 




m . + m . , 
J J-1 
2 


m . + m , . 

j j-1 
2 


(4.14b) 


As i® outside point for j = 1 and m^ is so for j 


NY. 


we use 


(4.14c) 


m. = 


^i.j±l/2 


n.j±i - n.j 


(4.14d) 


From equation (4.14b), (4.14c) and (4.14d) we get 


’^i.j-l i,j-l ’^x.j + l i,j + l 

Similarly, for = simulation depth 

n-M,j "i + l,j = ""i-l.j “i-l.j 


(4.14e) 


(4.14f) 


By equation (4.146) and (4.14f) we replace points outside 
domain in equation (4.11). This method is called image fcoint 
technique . 

(iii) Now comes the interface between Si and Si02, under 
gate. Special care is taken in discretizing the Poisson’s 
equation at interface. 

Boundary condition for interface is given by equation (2.16) 
of chapter 2. It ia rewritten for discretization here. 


^ = . f 

Oyf 


(4.15) 


^*Jsi 1 

Sx 




at i = 1 under gate. 



Expanding the point (i, j) about (i+1/2, j) and using 
equation (4.5), we get, 


/dy 


= - 



\ax 

dx/i 

,jjSi 

^*^Ji-H/2,j ^ 

1 g. + 

^sx^ i , 


[dx i+1/2, j 






ay- 


Si 


-h (n-p-C) 


(4.16a) 


Similarly expanding (i, j) about (i-1/2, j) 


= 

\ax sxj ^ J I 


+ £E11 +^i-ip^v' + a^cA\ 

idx dxl , ^ , 2 l"T 2^ 21 . 

i-1/2, j ^ax dv^Ji,j 


SiO, 


' J I "■*•''2 

(4.16b) 

The point (i-1/2, j) is inside the oxide where Laplace 
Equation is valid, 


i7J 


rl (V + «) = 0 


(4.16c) 


Since, oxide thickness (t^^) is very small, one dimensional 
linear potential drop inside the oxide is considered. So, the 1®^ 
term in the r.h.s. o£ equation (4.16b) is approximated by 


^ ^ ~ . j ~ ^i . J 

0x d>x t 

OX 


(4.16d) 



where 


^GCapplied) ^Flat Band 


(4.16e) 


and 


Flat Band 


= ^HS - 


'INT 


ox 


C4. 16f ) 


where, : Gate metal - Semiconductor work function difference 

: Gate oxide capacitance 


With this definition of v'q» we don’t need to treat Qjux 
equation (4.15) separately, we drop this term from equation 
(4.15) and proceed. Using equation (4.16d) and Laplace equation, 
we get. 




- n,j 


Oi 


i, j 


ox 


oxf 

■2"r 



^Ji, j|Si02 


(4.16g) 


Now, we put equation (4.16a) and equation (4.16g) in 

modified equation (4.15) is dropped with v'q definition], 

and we get after discretizing etc. terms as done earlier, the 

9yL 

difference equation at interface ; 




[at i = 1, j s under gate] (4.17) 


•4.3b Dlscratlzatlon ot Contlmiity Equation i- 

The major difficulty of the discretization of continuity 
equation is to find a proper numerically stable formulation using 
only the information at the mesh points. The naive discretization 
proposed by Van Dorpe [4] has been proved to be unstable [3], [8], 
[32] . It is also tested by this program, and is found to be very 
prone to overflow traps and inconsistent solution when large 
voltage is applied and mesh spacing is not very close. It is 
clear from difference scheme of Poisson’s equation that ^ is 
interpolated linearly between neighbouring mesh points. A naive 
linear interpolation (used in [4]) of n and p is surely 
inappropriate as n and p vary as the exponential of yf, unless the 
variation of y between mesh points is less 2KT/q [3]. 




So, we’ve chosen the discretization scheme proposed by 
Scharfetter and Gummel [3]. Although, they have proposed it for 
one-dimensional Read-diode oscillator, it is extended to 2-D as 
per the requirement of the simulation here. 

The outline for discretizing the continuity equation for 
electron is discussed here for hole it can be derived analogously. 

The normalized electron current density expression reads 


j - nC-^) + V Vn 

n *^n ^ ^n 


(4.18) 


Uhere y is normalized mobility, y = u fu 
is broken in two orthogonal components, 

J = J + J 
n nx ny 


(4.19a) 


where , 


J 

nx 


y n(- 
' n 


5x 


) + y 


Sa 
n ^x 


(4.19b) 


ny 


= y n(- 
n 




) + r 


dn 
n dy 


(4.19c) 


Now, we start our derivation with equation (4.19a) rewritten 

as , 



(4.20) 



(4.21) 


where, E 


Ue concentrate, on the path between two mesh points (i 


and (i+1, j) as shown in figure 4.3. 


(i, J) 


Ci+l/2,j) Ci+l,j) 


Figure 4.3 


U : a local variable at x = x^, U = 0 and at x = *i+l» ^“^i 

A ^ ^ 

" m 

If we consider E, r and J to be constant to the 

’ *^n n 

X 

mid-vertex (i+l/2,j) value over this path, equation (4.20) 
becomes a Leibnitz differential equation as. 


. E n = 

m * “1+1/2, j “ X 

"i+1/2, j 


(4.22) 


with boundary conditions. 


at U = 0 


n = n 


at U = hj^ , n = *^i+i/ 2 , j 


(4.23) 


Solution of equation (4.22) (without Boundary condition) 






.i = r “*i + l/2,j ^''*^1 + 1/2, j 

— P , , ■* # +1 5 


(4.24a) 


‘i + 1/2, j 



Where, C is a constant and evaluated by equation (4.23) and 


rearranging the final expression, we get. 


T _ B i+l,j 


n, • 
1 » J 


.-exp( + E^^^^2,j*^i> 


(4.24b) 


Equation (4.21) discretized as : 


1 + 1/2, j 


_ ^i,j - ^i+l,j 


(4.24c) 


Combining equation (4.24b) and equation (.4. 24c), we get. 


y 

'i + l/2,j" Tf^[*'i + 1, j®^^i + l, j " "^i, j ■ '^i + l,j^] 


(4.25) 


and J 


where, B(t) = — s : Bernouli function 

e"*^ - 1 

Similar will be the expressions for J J 

1-1/2, j ^yi,}+l/2 

^"^i, j-1/2 

Now, we consider the steady state continuity equation for 


electron. 


aj S3 

nx . ny _ B 

Sir- * Sir- - * 


(4.26) 



[In this simulation we neglect R, as the MOS not in 
breakdown, for the shake of completeness we keep R here]. 

Discretizing equation (4.26) as done in case of Poisson’s 
equation. 


nx 


i + 1/2, j 


h. + h. , 
1 1-1 


”"1-1/2 .j , 


OCh) 




ny 


i, j + 1/2 


^^i, j-1/2 


m . + m . . 
J J-1 

2 


0Cm)f-5-^^l =R. . 

Uy Ji^j i.j 


(4.27) 


Neglecting L.T.E. of order 0(h), 0(K) in equation (4.27) and 
combining equation (4.25) and similar expressions, we get. 
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1 + 1/2. j 
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1-1/2, j 

1, j+1/2 

1, j-1/2 

= R , 

* f J 


l^i-lClii + 
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mj(mj + mj_j^) 
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"i.j ■ ''l.J-1^ ' 




n«j_iC«j + «»j_i) 
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(4.28) 



where. 




+ y ) 

i + l,j *^1,3 


and L.T.E. in equation (4.28) is 


n 



aj 


dJ 

< 0(h) 

nx 

+ 0(m) 

ny 

dx 

Sy 


(4.29) 


(4.30) 


Analogously, difference equation for hole continuity 
equation can be derived and the final form of equation will be 
same as equation (4.28), except n, ,’s are replaced by p, ,’s, 
subscript of y is changed to p and the argument of Bernouli 
function are multiplied by -1. 

Now, boundary condition for continuity equation is to be 
imposed . 

(i) At the ohmic contacts, n and p are known from equation 
(2.17). 


(ii) At the fictitious boundaries we use image point 
technique to eliminate points outside simulation domain, as : 

”i,j-l ^ ”i,j + l ’ + l J ^ (4.31a) 

and ^i'l'l J ™ ^i""l J * ^1*^1 j ^i"*! j ^ *” (4.31b) 


(iii) At the Interface 



We impose no current component normal to the interface. 


equation (2.18) reads 


0 *>^<1 Jpjj, = 0 at i = 1 , j : under gate 


Without the loss of generality, we assume, h^^ = ^i-1’ 

h^_j^ is fictitious at i = 1 , and we write (for electron) 


J + J 

1+1/2, J 1-1/2, J 


(4.32) 


As J = 0 at i = 1 , J : under gate, equation (4.32) 

“*i, j 


implies , 


“'"i-l/2,j ”*i + l/2,j 

Combining equation (4.33), (4.27) and (4.25) and similar 


(4.33) 


expressions, we get 


'■”i + l,j ®^^i+l,j ~.’^i,j^ “ '^i.j ~ ^i+l,j^^ 


i+1/2, j 


‘i, j + 1/2 


^i.j + 1 ^^^i,j + l ~ ^i,j^ ~ ^i,j ~ n,j + l^ 

mj(m^ + ®»j_j^) 


‘i, j-1/2 


°i.i ~ n.j-P - - n.jJ 


(4.34) 



The difference equation for hole continuity equation at 
interface is analogous. 

4 ,. 4 , COMPUTER IMPLEMENTATION OF BERNOULI FUNCTION s- 


The Bernouli function used in discretizing continuity 
equations, defined as, 


BCt) = -_-l_ 
e - 1 


C4.35) 


In the actual programming of B(t) special care is taken to 
avoid floating point overflow and underflow traps. it is clear 
from equation (4.35), at t = 0 , B(0) gives divide by zero error at 
t » 1 and t « 0 exponential terms may give overflow and underflow 
error. So Bernouli function is implemented as suggested by Hart, 
Cheney, Lauson and Maehly [33], [25]. This implementation is 
basically defining B(t) by different expressions at different 



when t :£ tl 
when tl < t < t2 

when t2 < t < t3 

when t3 < t < t4 

when t4 £ t i t5 


0 


when t5 < t 


(4.36) 



The 


constants tl to t5 


depend on individual computer 


hardware. 


For tl, 


For t2. 


For t3. 


For t4, 


For t5. 


They are defined by 



1 




[t2 < 0] 



t3e 


-t3 


1 - e 


-t3' 


[t3 > 0] 


1 


-t4 


e 


* ^ f 


1 


-ts 


(4.37a) 

(4 .37b) 

(4.37c) 

(4.37d) 

(4.37e) 


The equal to (=) sign is quoted to mean that it is not 
mathematically equal, but the numerical values calculated by the 
computer are equal on both sides of the equations (3.38a...e). 

The constants evaluated for HP 9000 computer are as 


tl ac - 14.4954662 
t2 - 0.650278232 x 10“ 
t3 as + 0.6083747372 x 10 
t4 ac 14.4954662 


t5 a. 87.3300171 



4 ,. S ADAPTIVE GRID GENERATION : - 

As the discretization errors depend on the distribution of 
the quantities ¥*» n, p, a suitable mesh cannot be estimated a 
priori, i.e., without the knowledge of the solution. Therefore, 
grid generation is performed adaptively, a preliminary solution is 
calculated on the basis of an initial mesh (in sec. 4.2), then the 
mesh is adapted to this solution and again basic equation are 
solved. Re-generation of mesh is repeated whenever necessary. It 
is to be mentioned here that, the preliminary solution of basic 
equation need not be very accurate solution; it, even, may be the 
initial guess followed by solution of poisson’s equation and p & n 
adjustment . 

An easy method, used in adapting the mesh is described below 
From theory we know v' varies as an exponential function of 
spatial coordinate away from the interface of Si - Si02 when there 
is a gate voltage. The function is crudely of the form (only 
x-variation) 

V'(x) = A e (4.38) 

where A, B are constants. 

Expanding equation (4.38) about neighbouring mesh point, we 
get, after some rearrangement. 


h^B^ 


V>(x * h) - y(x) 

yiTxli 


hB + 


2 


(4.39) 



Ue use L.H.S. of equation (4.39) as maxiinum relative error 

n't 

The 1 term R.H.S. of (4.39) describes the variation as 
linear one between two mesh points; this assumption is made during 
discretization of Poisson’s equation (Sec. 4.3). So the local 


2 2 

error in linear interpolation of ^ is given by }h B /2j. 
^max 2 

Now (hB) is calculated from equation (4.40a) and 
equation (4 . 40b) [derived from equation (4.38)], 


(4.40a) 
feed to 


V'(x) - v^Cx + h)_ ^ .-hB 


(4.40b) 


Equation (4.40a) and equation (4.40b) gives an estimate of 

maximum allowed fractional change in normalized yf between mesh 

-3 

points, for which local error is less than , 

fractional change in normalized is calculated to be 4.39%. 
Quantitatively speaking, if we apply 6V at the maximum and built 
in potential is within 0.9V, then one grid is to be inserted 
whenever there is change of 11.6 thermo-volt between mesh points. 

In discretization of continuity equations (Sec. 4.3b), we 
considered the variation of V' linear between mesh points, so this 
rule of grid adaptation is valid for n and p distribution also. 

The initial coarse grid is generated accounting for the bias 
values, doping profile and important boundaries. 

Ih x-direction . Si - Si02 boundary is very important, as 
inversion layer (if formed) is limited within few thermo - voltage 



& *fc 

fall from interface. So, I used 1 grid spacing as the depth 
within which yf • drops by one thermo-voltage. Then spacing 
increases exponentially as we move away from surface. 

In y-direct ion . fictitious vertical boundaries are prone to 
large local truncation errors, so y-grids are closely spaced 
there. To track the lateral diffusion of source and drain, grids 
are inserted for every 100/CC fall in doping concentration. this 
thing is done in x-direction upto the junction depths. The 
initial grid spacing is sinusoidally varying in the effective 
channel . 

In this way minimum number of grids (18 x 20) for a 
particular case, and maximum allowable grid adaptation is done 
upto (49 X 49). 





ALGORITHMS TO EXTRACT RESULTS 

In the preceding chapter, three coupled a emi conductor 
equations are linearized as required and discretized. A large 
system of linear simultaneous algebraic equations, with the values 


of the dependent 

variables of 

the differential 

equations 

at 

discrete points as 

unknowns, are 

obtained . 

Direct 

solution 

of 

this system becomes 

impossible . 

For example, 

if grid 

points 

are 


Nx.Ny = 20 X 20, the number of unknowns becomes {3(Nx.Ny)} and the 
size of the coefficient matrix will be 3(Nx.Ny) times 3CNx.Ny). 
Uhich amounts to (1200, 1200) at the minimum number grids. 

To solve this large system, GUMMEL’S ALGORITHM [1] is 
popularly used in device simulation. 

5.1 GUMMEL'S ALGORITHM s- 
5.1. a The Basic Idea : 

Gummel’s idea was, in essence, to solve the coupled 
semiconductor equations by independently solving each equation, 
assuming the knowledge of the result of the other equations. The 
steps involved are as follows : 

(i) Prepare initial guess for y. n P- 

(ii) Assuming knowledge of v to be exact, solve discretized 
electron continuity equation to get better approximation to n 

(iii ) Assuming knowledge ^ to be exact, solve discretized hole 
continuity equation to get better approximation to p. 



(iv) Now, assuming knowledge o£ p and n from steps (ii) and (iii) 
to be exact, Poisson’s equation is solved, to get better 
approximation to y- 

(v) Test the convergence and eventually return to step (ii). 

The effort to solve for one cycle of this iteration, is 

obviously less than any other solution methods like Newton’s 
method, because the rank of the matrix for one decoupled equation 
is one-third of the rank of the total system. 

5. 1 • b Convergerice i 

A complete theoretical proof for the convergence of the 
algorithm has not been published yet. However, strong theoretical 
indication has been given by Hock [35]. The ’’practical” value of 
this algorithm is not underestimated, many authors [Vandorpe [4], 
Scharfetter and Gummel [3], Mock [6], Selberherr et al . [13]# De 

la Moneda [19] and several others] have used Gummel ’s algorithm 
with excellent success as they claimed. 

6. 1 . c Modification i 

Gummels algorithm is slightly modified and simplified for 
level 2 model [Sec. 3.1], used in this simulation. In level 2, 
only majority carrier continuity equation is solved. So, any one 
of the steps, (ii) or (iii), is replaced by corresponding minority 


carrier adjustment, as: 



(5.1a) 


n 


K + 1 _ K+1. 

- n expCoi ) 


Or ..K + 1 K , K+1., 

• P = P expC-a ) (5.1b) 

superacript k is the value at K-th iteration, and a is 
®rror in ^ as described ago, and y is adjusted as 
. K + 1 


K ^ K 
^ + a 


(5.1c) 


Initial Guasa t 


In strong inversion region it is reported, the decoupled 
needs 200 or more iteration [16], [19]. One reason for 
^^ia, is due to the strong coupling between continuity and 
^Piason's equations, under strong inversion. Other reason is the 
^a.rg« error in "thte initial guess. It has also been tested in this 
®i*»Ulatlon program that crude initial guess, sometimes caused 
‘divergence and floating point overflow. So, extensive care is 
i^aken in the preparation of the initial guess. For, the 
'Uniqueness of the solution [Mock [34]], specified initial data (y/ 
*'* P) which fulfill the Poisson’s equation is necessary. 

ui i; 1 shows the flow chart of initial guess preparation. 

Increment in the terminal bias is set 20V,j. by trial and error 
®®tho<t to avoid over flow in equation (5.1a) and equation (5.1b), 
increment is greater than final bias value, directly final bias 

Set. 

Figure 5 2 in flow chart for modified Gummel’s algorithm, for 
^ Single operating point. For calculating I-V characteristics 
thia ia repeated for operating each point. 



Processed input data on the 
generated grid points 


Set increment in terminal biases. 
Set IMREF’s in each impurity region 
equal to corresponding bias value 


Set n, p to thermal equilibrium 

value and set y as : 

^-An <pp~v 

n = e ^ , p = e^ ^ 



Increment terminal biases 


/tinal\ 

value 

reached 


To original Gummel’s loop 


Figure 5.1 : Flow Chart for Initial Guess 














Input guess v'» ti. P from 
initial guess routine 



Figure 5.2 : Flow Chart for Gummel’a Algorithm 












S. a SOLUTIOM OF LARGE SPARSE LINEAR SYSTEMS 


The discretized (5-point molecule) equations 
and continuity equations are rearranged 
implementation, with coefficients in the form of 
spacings, which is necessary to minimize the 
truncation error during matrix inversion. They 
following 


for Poisson’s 
for computer 
ratio of grid 
floating point 
look like the 


Poisson’s equation 


"i,j-l^ mJZl ^ “i-1 


m.+m._. f . 

,j^ ^ h^ ^ + m._^)C^ + ) 


1 i “s-i ) 
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. j j“ls . ^ i i-1. 
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j-1^ m 




i + l,j' h 




(5. 2a) 


Electron Continuity equation 
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(5.2b) 


Hole continuity equation will look a like so, it is not written 
here. 

The family of equation (5.2) gives rise to a large sparse 
system of linear system as : 


(MX.NY)(NX.NY) ^ (NX .NY) ( 1 ) ~ ^ ^ (NX.NY) . ( 1 ) 

The coefficient Matrix has at most 5 non-zero elements. The 
direct solution of this system is very costly in terms of computer 
memory. Iterative schemes are found the most suitable for 
solution of these semiconductor problems. 

The main disadvantage of iterative scheme is their alow 



convergence. A relative merits and results of different iterative 
schemes [ point- Jacobi , S.O.R., S.L.O.R., ADI etc.] have been 
discussed in [42]. 

As Poisson's equation is solved quite a number of times in 
the algorithm than the continuity equation, so a fast algorithm 
is chosen. Stone’s strongly implicit procedure (S-I-P). 

S. 2. a Stona's S-I-P s 

Stone’s idea [36] was to perturb the original [A] by another 
matrix [N] of same rank but the norm of [N] is much smaller than 
that of [A], so as to factorize [A+N] easily and factorized lower 
[L] and upper [U] triangular matrix gives only non-zero elements 
corresponding to the non-zero places of the original [A] . Figure 
5.3 shows the non-zero diagonals of [A], [L] and [U]. 

In figure 5.4, solid circles represent the standard 5 point 
molecule for discretization scheme. To get the effect of [N] , 
perturbation matrix, we include the effect of crossed (X) points. 
Expanding the (X) points in Taylor series and neglecting higher 
order terms : 


T a -T 

^i-i,j+i n,j 

+ T 

j+l 

+ T 

i-l, j 


(5.4a) 

^1+1, j-1 i,j 

T 

^i+i,j 

+ T 


(5. 4b) 

To ensure the effect of 

perturbation 

minimum 

right 

hand 

sides of equation 5.4 

are multiplied by oi (0 

< « < 1) 

and 

added 


to left hand aide and finally incorporated in the left hand side 
of discrete equation [A1[T] = [q], L.H.S. becomes 
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••trice* for Stone's S-I-P 



Fig. 5.4 

Discretization eolecule for the perturbation of 
^^^ffient ■atrix CA3 in Stone's 8-I-P 
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B 


ij 


* '■i, j[’'i + l. j-l ■ ■*■ ’^i + l.j * 


T 4-H T 

C5.5) 


Coefficient ^ ^ through ^ are as shown in figure 5.3. 

The coefficient matrix now becomes 

[A' ] = [A + N] = [L3[U] (5.6) 

Collecting the terms in equation 5.5 and comparing with the 

R.H-S. o£ equation 5,6, we get 

^i,j ' ‘'i.j ®i,j-i 

®i,j ^i.j ^i-i,j 

The terms multiplying them are partially canceled dependent 

variable at (i-1, j+1) and (i+l, j-1). 

The relation obtained relating coefficient of [A+N] and 

[L] X tU3 ate 




(5.7a) 


*>i.j “ ®i,j -“S,j 
"i.i “ “i,j ' “°i.j 
‘"i.j * ”1,3 “i-i.j 


I . . + aC, . + «G. , 

i. J i. J i, j 


‘^i.j ®i,j “ ^i,j 



"^i.j ^i,j " “i,j 

■ “°i,J 

(5.7b) 


where, the lower case letters denote the coefficients of [L3 and 
CU3 and capital letters for [A+N3 , as shown in figure 5.3. 

From equation 5.3 we can write 



[A+N][T] = [A+N][T] - C[A][T] - [q]) 

It gives the iterative scheme equation 5 . Sb 

= [A+H][T]‘= - ([AKT]*^ - Cq)) 


At any iteration step, coefficient of [L] and [U], i.e., b, 
c» d,, e, f are calculated using the equation family (5.7) by 
elimination method and equation (5.8b) is solved by LU 
factorization, until convergence occurs. 

The parameter «, called iteration parameter, is chosen 
optimally. The range of « is [0, 1]. But two small value of « 
causes convergence problem. A set of a. values ranging from 0.6 
to 1.0 is tried; a value around a. - 0.98 is found suitable for 
faster convergence of S-I-P. 

S. b S^c«rssiv« Line Relaxation s 


Continuity equation yields a coefficient matrix which is 
rather ill conditioned. It has been found very difficult to solve 
continuity equation by S-I-P. So, this successive line relaxation 
Iterative technique is used here. 

It is baaed on a sequential update of the solution according 
to the block of nodes [37]. For these techniques unknown vector 
[T] is divided into NX blocks corresponding to the horizontal 
lines of a grid structure (figure 4.1). The elements of 
CAJ,CT3.Cq] becomes blocks of matrices as [B]^, [C]^, [D]^, [T]^, 
and {q3^ for the i^*^ X-line on the grid (figure 4.1), and [A] 



becomes block tridiagonal matrix [37], as 




C 

[A.] - n 


The diagonal blocks are again tri-diagonal and [B]^ and 

are only diagonal matrices. 

It is solve by the iterative scheme as 

=Cl-w)[T]^+u[C]]^^|[q]^ - [B]. [T]^_^ - [D]. (5.9) 

for i = 1 to NX 

Where, the relaxation parameter, w, is in the range 0<w<2 
[371, [381. 

As done in case of a, value of w is set to 1.2. 

This type of iterative scheme requires low amount memory and 
they are very easy to program. 

S. 3 pcsrr “PROCESSING C«F V, n, p i- 

Ue get n, p as solution of semiconductor equation by 
Guamel’s algorithm; now we extract from the results 

obtained . 

B.3,a Current Calciilation i 

Pao and Sah's model [40) was developed long ago to calculate 




current in MOS, excluding the effect of source and drain 
junctions and taking one dimensional Poisson’s equation in the 
channel. This model is modified by De la Moneda [19], is used 
here . 

Assuming predominant current flow in lateral direction, for 
electron, it is re-written as 


ny 


y e 
^ n 


y/ de 


•<Pn 


^*5 


(5.10) 


Neglecting recombination generation current, and integrating 
the differential element [figure 5.5] dx dj^ over X-Z plane, 
we get 



where d and U are defined in the figure 5.5. 


d 

Defining F (■ 5 ) - J ''n 

o 


(5.11b) 


and carrying out an additional integration over channel length L 


^-<^nCo) ^ ^-^n(L) (5.11c) 

T » t| 1 

0 

Where x « 0 : source contact edge near gate 
X * L : drain contact edge near gate 
Siitilar Is the current relation for p-channel MOS. 
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5. 3» h Threshold Voltage CV ]> 

th 

1 « 




^th’ simulation and 

practical 

problems , 

is 

the gate 

voltage at which the device 

sinks 0.1 

pA X (W/L) 

amount of 

current [41]. Usually V . is 

th 

calculated 

f or V = 

DS 

0.2 

V or any 

value close to it. at higher 

V V 

DS' ''th 

decreases 

due 

to drain 


induced barrier lowering [18], 


One scheme, derived from the subthreshold current expression 
of nos, which varies exponentially with V__, states [14] 


th 


rJ-l 

th 


+ In 


f -^1 


1^1 

{ j 


-1 


(5.12) 


where, superscripts are the values at j-th iteration of the 
Newton-Rapson scheme given in equation 5.12. 


I , = 0.1 (U/L) M 

th 

'^th * starting guess for ^^th from long channel 

'^th *^®*^®**^* obtained in any text [10]) 
drain current at Vg « 

The derivative La calculated by giving increment on 
+ AV (fti 50 mV). 

Cl til 

The stopping criterion for convergence is set 


S 1 mV. 



5.4, CONVERGENCE CRITERION OF GUNNEL* E ALGORITHN USED *- 

Figure 5.6 shows the convergence criterion of the iterative 
algorithm. The normalized error is the root mean square error in 
potential at i-th iteration. Iteration counts are the Iteration 
number of the Gummel’s algorithm and counting starts after the 
initial guess loop ends [figure 5.1]. It is found that 10 to 35 
iterations through Poisson’s equation solution and p“U adjustment 
and bias increments are needed for a good initial guess; the 
number depends on bias values and number of grids. 

It is found in figure 5.6, when the MOS not in strong 
inversion (Vg = 0.5V, = 3V) only 7 Gummel’s iterations are 
needed for convergence. But under strong inversion (Vg = 3V and 
Vjj = 3V:same) convergence is slower, takes 60 iterations. This is 
due to strong coupling between Poisson’s and continuity equations. 

It is also noted that, in this simulation, initial guess 
iteration plus Gummel’s iteration is within 100, owing to the 


better initial guess in this simulation. 



Convergence criteria 
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SIMULATED EXAMPLES 

It is rather difficult to present all the possible simulation 
results here. Examples are chosen which cover some of the problems 
and remedies associated with the miniaturized MOSFETs . The 
examples are simulated with the help of the two dimensional 
simulation program developed . 


6.1 SIMULATED DEVICES i- 

In this series of examples the devices chosen are the short 
channel MOSFETs with no implantation in the channel ( DEVICE 1), 
with one shallow channel implantation ( DEVICE 2) and with two 
channel implantations ( DEVICE 3) . 

6.1a. DEVICE 1 i 


Figure 6.1a shows the input file [ explained in the appendix] 
of the DEVICE 1, n-channel MOSFET. The gate mask length (1) is 1.5 
)jm, width (w) is 10 fjm, oxide thickness (tox) of 500 A, source and 
drain contact lengths (upto which simulation is done), are 0.3 

15 

Source and drain are implanted with a Phosphorus- dose of 10 
2 

atoms/cm , range (rp) and straggle (drp) are set so that junction 
depth be approximately 0.3/jm; the substrate concentration (subs) 
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title DEVICE-1 

device type«'n', *d=0.3e-4, l“1.50e-4, w«10.e-4, 

+ tox*500.e-8, d*£.5e-4 

bias ud»3.0, ug®0.0, ub=0.0, ub*0.0, vfb=-0.8 

comnt No channel implantation 

profile subs*-l.el5, dose® 1 . OOel 5 , rp®0.0, drp*6.01e-6 
comnt Distribution of potential, n, p 4 doping profile 
option model='distrib' 
end 


title DEVICE-E 

device type*'n', sd®0.3e-4, l®1.50e-4, w»10.e-4, 

+ tox®500.e-8, d®£.5e-4 

bias ud»3.0, ug=0.0, us*0.0, ub=0.0, vfb=-0.8 
comnt one shallow channel implantation 
profile subs=-1.e15, dose® 1 . OOel 5, rp=0.0, drp®6.01e-6 
+ dosed =-3 . el 1 , rpc1=0.0, drpc1*6.e-6 

comnt Distribution of potential, n, p 4 doping profile 
option model* ' dist rib ' 
end 


title 

device 

+ 

bias 

comnt 

profile 

+ 

+ 

comnt 
opt ion 
end 


DEVICE-3 

type*'n', sd*0.3e-4, l®1.50e-4, w®10.e-4, 

tox*500.e-8, d=£.5e-4 
ud=3.0, ug®0.0, us*0.0, ub®0.0, vfb®-0.8 
two channel implantations 

subs®-1.el5, dose=1 . OOelS, rp®0.0, drp®6.01e-6 
dosed ®-3. et t , rpd®0.0, drpd=6.e-6 
dosec£*-£ . el 1 , rpc£*3.9e-5, drpc£®8.e-6 
Distribution of potential, n, p 4 doping profile 
model®' di St rib ' 


Fig. 6.1 

Input files for <a) DEVICE 1, (b) DEVICE £ 4 

(c) DEVICE 3 
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15 

is 10 /cc of p-type impurity (Boron) .The device is biased at : 

us=ub=0 V, ue= Ov, ud= 3v, flatband voltage is set to -O.Sv (only 

taking work function difference of Si and gate metal Al). 

From the quasi 3-dimensional plot of the doping profile set 

by the program (fig. 6. 2a, source on the left) source - drain peak 

20 

concentration is measured to be slightly less than 10 /cc (value 
given in the output file is 0.6655B20 /cc).The output file records 
vertical junction depth (Xj) equal to 0.2d/jm and lateral spread is 
approximately the same. So the metallurgical channel length is 
reduced to 0.94pm . 

Figure 6.2b shows the quasi three dimensional and figure 6.2c 
shows the contour plots of the potential distribution within this 
DEVICE 1 . The potential at source and drain contacts are found to 
be higher than what is applied by O.Sv, which is built in 
potential of the junctions. In the depletion region of the reverse 
biased drain-bulk junction, it decreases monotonically , and it is 
nearly constant in the highly doped source & drain regions . It is 
clear from the source-channel potential variation that there is 
hardly any barrier. So this device is on even at zero gate 
voltage . 

The above feature is very clear from the electron 

distribution ( fig. 6. 2d) . The surface concentration of electrons 

15 

in the channel is high enough (= 6x10 /cc, above the substrate 
concentration ) to cause inversion in the channel . It is also to 
be noted that surface concentration of electron decreases near 
drain end; it explains the occurrence of pinch-off in the channel. 
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One more thing very clear from fig 6. 2d, is that, there is an 
unexpected bulge of electron concentration at a depth of about 0.5 
to 0.6 ^m . This is the proof of ’’punch-through” present at the 
said drain voltage, the bulge of electron supports the 
punch-through current flowing away from the surface. The 
punch-through is further supported by the merge of drain and 
source depletion regions in the figure 6.2c. This punch-through 
and unexpected formation of inversion layer (i.e., decrease of 
threshold voltage ) are the major problems of short channel MOS. 


6.1b. DEVICE 2 i 

The input file of this device is shown in figure 6.1b. All 
other conditions are same as DEVICE 1 except, there is a shallow 
p-type channel implant to increase the threshold voltage, i.e., to 
make DEVICE 2 normally off. 

From the impurity distribution (fig. 6.3a) peak concentration 

1 6 

and the depth of the channel implant are measured to be 2x10 /cc 
and 0.257^m respectively. This high concentration inhibits 
formation of inversion layer. 

Figure 6.3b and 6.3c show the quasi 3-d and contour plots of 
the potential distribution in DEVICE 2, under same bias as before. 
The quasi 3-d plot looks like that of DEVICE 1, but the contour 
plots show there is slight barrier formation in source - channel 
compared to DEVICE 1. It also shows that the single channel 
implantation is unable to separate the depletion regions of the 







esiOiio *0 ^ V 
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source and the drain. So punch- through occurs here. This is 

further supported by the similar bulge of electron concentration 

away from the surface (fig. 6.3d) . 

The surface concentration of electrons has decreased to the 
14 

value 1.2x10 /cc ( less than surface accepter concentration) as 
expected by the channel implant, no inversion occurs at ug=0v. So 
one of the short channel problem has vanished, but there is still 
the problem of punch-through. 


Figure 

6 . 3e is 

the quasi 3-d 

plot of 

electrons 

looking 

through the 

gate . 

It exaggerates the 

pinch-off 

of the 

channel . 

Near drain. 

the surface electron concentration 

sharply 

decreases 

at the surface but 

it increases at 

a slight 

depth 

from the 

surface. So, 

in the 

pinch - off region 

, current 

flows away from 

the surface 

and comes back to the drain contact 

. Hence 

we need 


2-dimensional simulation to track this phenomenon . 

6.1c. DEVICE 3 s 

As shown in figure 6.1c, DEVICE 3 has one more p-type deep 

implantation, which we will see controls punch-through, while 

other input conditions are same as that of DEVICE 2. The 

concentration profile ( fig. 6.4a ) gives the peak concentration 

1 6 

of the second implant is 1.07x10 /cc occurs near a depth of 0.4^m 
and spreads upto 0 . 6jum inside the bulk. It can be seen from the 
figure 6.4b that there is a source channel potential barrier at 
the source end. Two distinct depletion regions, for source and 
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drain, can be noted in the isopotential plot (fig. g 

explains the clear suppression of the punch-through . 

Figure 6.3d is the electron distribution in DEVIce 3 The 

major difference from the last two devices, is the absa^ 

“®bce of any 

bulge of electron concentration away from the surface 

this is one 

more indication for the suppression of punch-through, ti... , 

^he surface 

concentration of electron is less than ®*^t'fac6 dopant 

concentration, so no formation of inversion layer at 

zero gate 

voltage . Hence we can infer that, two short channel effects low 
threshold voltage and punch-through, are overcome in DEVICE 3 

So far, no distribution is shown for HOSFET strong 

inversion. Figure 6.4e and figure 6.4f show the potential and the 
electron concentration distribution under strong i^^versiou ( 2v at 
gate ) respectively. The concentration of electron at the surface 
is very high ( clear from comparing with fig. 6 . 4d ); decrease 

of electron concentration near drain contact shows pinch off as ud 
is set at 5v. 


6.2 . COMPARISON OF TERMINAL CHARACTERISTICS OF THE THREE DEVICES 


DESCRIBED t- 


So far comparison is made on basis of the it>ut ion of 
potential and majority carrier concentration of DEVICE DEVICE 
2and DEVICE 3 at a particular operating point . Nov comparison 
will be shown with respect to the simulate^ terminal 


characteristics for several operating points . 
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6. Sa. SSub'threshold charac'tarls'blcs : 

It la the drain current va . gate voltage characteriatics 

below the threahold voltage. The subthreahold region is 

particularly important in low power and low voltage applications, 

such aa , MOS used as a switch in the digital logic circuits. These 

characteristics say how the MOS -switch is turned on and off. 

Figure 6.5a shows the subthreshold characteristics for the 

three devices at 0.2v of drain bias ( 0 . 2v is chosen as low state 

logic voltage is around this value). It is clear from the semilog 

plot that, the drain currents fall off exponentially with gate 

voltage, which tally with the theory [10]. The slight shift from 

the linear fall of the semilog plot at very low current below 
-13 

10 A is due to the accuracy limitation of the simulation 
program. 

One important parameter extracted from the plot is the 

subt?xr'&s?u>ldt sioin^ (S), defined conventionally as the inverse of 

the elope of (in los^^ ecale) va. gate voltage. S aeaaurea how 

fast the MOS - switch can be turned off and on, the steeper the 

slope (hence smaller S ), the less gate voltage sweep is needed to 

change state of the switch. For V^g=0.2v, S is nearly same for all 

the three devices, approximately 80-85 mV/decade of 1^^ ; but the 

residual I. increases from DEVICE 3 to DEVICE 1. The drain 
ds 

currents levels off as gate voltage approaches threshold voltage. 

When drain voltage increases to 2v (fig. 6.5b ), there is no 



DRAIN INDUCED 

barrier lowering 



. , Fig. 6.6 

Surface potential distribution aionn 
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substantial change of S in case of DEVICE 3, but S increases to 
lOOinV/decade of in case of DEVICE 2 and a drastic change is 
observed in case of DEVICE 1(3 =200inV/decade ). The residual 
current increases the maximum in case of DEVICE 1. 

The shift of subthreshold characteristics at higher drain 
voltage is due to the dii^ain. truixkcefd. hax^x-i&x' (DIBL) effect 
[18] for short channel MOSFETs .This is explained in figure 6.6. 
For long channel Cl=!5.5/Lim) and for short channel Cl = 1.5pni) n-MOS, 
the surface potentials near source junction are shown. It is seen 
that potential barrier between source and channel decreases in 
case of short channel one ( both are biased at ug=0 & ud=5v ) . The 
reason for this barrier lowering is due the proximity of drain and 
source depletion region in case of short channel MOS. The leakage 
I^g is exponentially dependent on this barrier. 

6. 2b. Output characteristics s 


The output characteristics of the three devices for a gate 

bias of Iv are shown in figure 6.7. I. level decreases from 

ds 

DEVICE 1 to DEVICE 3, because surface electrons (inversion layer) 


are maximum in the 

1®* device and 

decreases 

by 

the 

channel 

implants in other 

two devices 

The slope 

of 

the 

output 


characteristics in the saturation region for the DEVICE 1 is large 
(so X is large) compared to other two devices; this is one more 
proof for the pronounced punch-through in DEVICE 1 . 



put 



Punchthrough characteristics 
Vgs=0v 

DEVICE 1 . 



Punch-through charactorist ic» of th* thro* d*vi 
Cug*0v3 




88 


6. ae. Pimeh-t.hrough characterlst-ics : 

Figure 6.8 shows the semi log plot of the punch-through 
characteristics of the three devices when gate voltage is zero. 
This plot shows the drastic increase of I , for V. crossing a 

CiS ds 

value, due to the flow of current along the merged depletion 
regions of source and drain. From these characteristics, 
punch-through voltage can be measured. is defined as the 

drain voltage at which channel current is equal to a reference 

~g — 9 

value. Different reference values in the range of 10 to 10 A 

— g 

are chosen for different applications [18] . Let us choose 10 A 
for the channel width of 10pm . DEVICE 1 is already on even at 
sero gate voltage, so cannot be measured here. For DEVICE 2 it 

Is found to be 1.4V . This effect of punch-through at low V^^ has 
>een observed in sec. 6. lb. 

For DEVICE 3, V is found to be 7.7V. It again proves the 
if f ectiveness of the deep channel implant in short channel 
ievices . 

£ec. Thr-AsHold volLaga a 

For the calculation of V^j^ by threshold model, we use 0.2V at 
train. Threshold voltage vs. channel length (1) is shown in 
•ig.6.9. It is found that typical short channel effect (fall of 
^th ^ starts at a channel length of 2pm approximately. 

DEVICE 1 is found normally on as its V^j^ is -0.006V for 



THRESHOLD VOLTAGE 



Fig. i.g 

Threshold voltage v« . channel length of the 




DEPLETION nMOS 






• ^ 05 * 


. Fig. 6.10a 

DapUtion n-HOSFET 

tin daplation mode 3 


DEPLETION nMOS 






Fig 6. 1 Ob 

Concantrat ion of alactron of 
HOSFET tin daplction aoda 3 


Oaplat ion 
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l = 1.5/Lini. For DEVICE 2 & DEVICE 3, they are 0.3V & 0.5V 

respectively. The effect of channel implantations on threshold 

voltage shift is pretty clear from this plot. 

6.3. SIMULATION OF A DEPLETION N-MOSFET 

A depletion type n-MOSFET of l = 1.5/jm, w=3^<m, tox=500A; the 

source-drain and deeper channel implant is similar to that of the 

DEVICE 3; only the shallow channel implant is changed to donor 

12 2 

type (dose = 1x10 /cm , rp=.01/Liro A drp=.06/jm), and the 

corresponding flatband voltage is -0.5V. 

This device is simulated in the depletion mode, with -IV gate 
voltage. Quasi three-dimensional plots of potential and electron 
distribution are shown in fig. 6. 10a and fig. 6. 10b respectively. 

The output characteristics is shown in fig. 6. 11 for gate 

voltage ranging from depletion to enhancement mode operation. The 
large slope of the characteristics is noticed, as there is no 
junction in drain & source to channel doping profile and the 
effect of drain voltage reaches close towards source through this 
continuous channel. 

Uith the examples presented so far, it is clear that the 
power of device simulation in understanding the behaviour and in 
the design of MOSFETs is enormous. 
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CONCLUSION 

A two dimensional numerical model for the MOSFETs has been 
described and results have been presented as extracted from the 
program developed for two dimensional planar MOSFET analysis. The 
results show that the device simulation program can simulate 
measurable quantities like output characteristics, threshold 
voltage etc, as well as the quantities which cannot be measured 
directly, such as potential distribution, carrier concentration 
distribution etc. within the device. These non-measurabl e 
quantities are important to understand the complex behaviour of 
the integrated device directly, useful in the design process of 
ICs. 


7.1. PROPCKSAL FOR FURTHER WORK s- 

Now, some scope for the modification and extension of the 
device simulator developed are outlined below . 

7.1a. ModlTlca-blon of* irhlB davlce simulator i 

There la, however, a need to develop more accurate model of 
the physical parameters, of the basic semiconductor 
equations, specially the doping profile. This may involve complete 
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re-evaluatioix of some commonly accepted assumptions and 
approximations . 

In the numerical model (level 2) used here does not include 
recombination-eeneratlon (R). Provision is kept for its inclusion 
with proper modelling, for the simulation of breakdown behaviour 
of MOSFETs. In that case both the current continuity equations are 
to be solved. So computation slows down, hence one can try with 
other algorithms (like s m&thodL ) to solve the non-linear 
coupled P.D.Es. 

One more modern trend is three dimensional device simulation. 
As the device dimensions are shrunk, the narrow width effect and 
effect of the neighbouring devices become predominant. Only the 
three dimensional simulations can characterize these effects. Then 
the memory requirement and the computing costs become the major 
problems; the main routines may be written in the assembly 
language for faster processing of data as done in case of MINIMOS 
[13]. 

7.1b. Integration of process, device & circuit simulators : 

In order to make any model predict the performance of the ICs 
prior to their actual fabrication, it is essential to integrate 
the process, device and circuit simulators. The device simulator 
cannot stand alone in the design and development of VLSI. The 
process parameters (like ion dose, implantation energy, oxidation 
temperature, etching environments etc.) are set for initial 
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design. These parameters are fed to the process simulators like 
SUPREM [28], SUPRA [29]. They give the doping profile for the 
device either in one dimension or two dimension. Then this doping 
profile is fed to the device simulator. The main problems to 
bridge the device and process simulator are of two types: 

If SUPREM is used, it gives one dimensional profile, but 
the device simulators require two dimensional profile. One has to 
pay attention to the extrapolation of 1-d profile to 2-d [32]. 

Ctiy The second problem is that the mesh of the process 
simulator is different from that of the device simulator. Error in 
the interpolation of the impurity concentration on device 
simulator mesh from the data on the process simulator mesh may 
cause double junction formation, which creates problem in device 
simulation [43]. 

Once these two problems are overcome, the gap between process 
and device simulator is bridged. 

The next is to integrate device and circuit simulations. 
SPICE circuit simulator uses analytical models[44] to describe I-V 
characteristics. The key to linking technology simulation to 
circuit is the ability to extract the model parameters from the 
I-V characteristics obtained by device simulation. This parameter 
extraction is a separate field. Several parameter extraction 
algorithms and softwares are developed so far [45]. These 
parameters are fed to the circuit simulator and small and large 
signal behaviour, transient response, noise performance etc are 


simulated. 
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For a full design of IC, several loops (depends on the 
complexity of the circuit ) are needed from process to circuit 
simulation and adjustment of the process parameters to achieve the 
final goal of the circuit . 

This type of integration can be done using this device 
simulator and two other softwares available at I IT, Knapur, SUFREM 
& SPICE . Only the bridging softwares are to be developed . Then 
we will be having a complete CAD tool for VLSI device design. 
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TITLE & END : 

DEVICE . PROFILE . BIAS . OPTION & COMNT 

The individual indents and the keys under them are described 
below. 

Ci!> DEVICE : It describes the channel type and the geometry 
of the planar MOSFET. The keys under this are: 

TYPE : for nMOS (pMOS) it is equated to ’n’ C ’p’)- 

L : gate mask length, W : channel width, TOX s oxide 
thickness, SD : source & drain contact lengths, D : simulation 
depth. They all are real numbers in cm units (need not be 
specified). All the keys are to be defined. 

Cit:> PROFILE : It describes the doping profile of MOSFET, 
from the parameters given here the program generates a 2d doping 
profile. The keys under this indent are : 

SUBS : substrate doping concentration (in /cc), if the 
impurity is p-type then the real value assigned to it must be 
negative (as ionized accepters are negatively charged ), for donor 
type it is positive. 

2 

: source & drain implant dose (in atoms/cm ), for accept 
or type value must be negative and vice versa. 

RP : range of source & drain implantation, DRP : straggle for 
the same. They are real numbers in cms. 

DOSECl : dose for 1st channel implantation (unit and sign 
convention is as DOSE), RPCl & DRPCl : are the range and straggle 
(in cm) for 1st. channel implantation. 
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DCkSECc! , RPC2 & DRPC2 : are the corresponding quantities for 
2nd. channel implantation. Any one or both of the sets (DOSECl, 
RPCl, DRPCl) and/or (D0SEC2, RPC2,DRPC2) may be omitted, then it 
is assumed that the corresponding channel implantation is not 
done. But the other four PROFILE keys must be defined . 

any BIAS : It describes the terminal biases and interface 
quality. The keys are : 

US, UG, UD & UB : are source, gate, drain and bulk voltages 
in volts respectively. 

PHIMS : gate metal-semiconductor work function difference in 

2 

volts, (JINT : interface trapped charge density (in coulomb/cm ), 
VFB : flat band voltage in volts. 

The value assignment to these keys are optional and if not 
specified default values are 0. If PHIMS & QINT are specified, no 
need to specify VFB. If VFB is specified, it overrides PHIIK; & 
CffNT. 


Civy OPTION : It describes what type of simulation the user 
wants and the parameters related to these simulations. The keys 
are : 

MODEL : it is equated to a particular option of simulation 
supported by the program. The options supported are as follows - 

’DISTRIB’ - gives distribution of potential, carriers & 
doping concentrations at grid points in the output file. 

’IVCHAR’ - gives the output characteristics ( I^^ vs for 
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different V ), two parameters are to defined along with this 

ga 

model. They are NGSTEP & NDSTEP , equated to the integers i.e., the 
number of gate voltage steps and drain voltage steps for output 
characteristics. The values assigned to UG & UD in the BIAS card 
become the final values for the output characteristics 
calculation. If MGSTEP & NEOTEP not specified, default is taken to 
be (1,1). 

’VTHRES’ - it calculates threshold voltage. UD of BIAS card 
must be a non-zero quantity and UG defines the initial guess for 
threshold voltage using the scheme described in sec. 5. 3b. 

’SUBTHRES’ - gives subthreshold characteristics. NGSTEP must 
be specified in this card, defines the number of points to be 
calculated upto the gate voltage specified by UG of BIAS card, and 
UD must be non-zero quantity. 

These model names must be quoted inside single inverted 
commas, written either in lower case or upper case letters, mixing 
of case is not allowed. 

T : this key is optional for specifying absolute temperature 
in “K, default value is 300 “k. 

COUNT : the card under this indent contains comment, if 
the user wants to give any, this is ignored by the program. 

Important points to note 

1) Uhen a line becomes too long the user may break the line 
into several lines and the continuation lines must start with a 
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plus sign (+) within first 10 columns of a line. 

2.) Except TITLE & END cards, other cards may be shuffled. 

3) All the lengths are in cms, concentration in /cc, biases 

o 2 

are in volts, temperature in K, doses in atoms/cm etc. as 
described above. One should not specify the units in the input 
file. 

4) A11 values are real numbers (written in F or E format ) 

except NGSTEP & NI3STEP, they are integers. Character strings 
equated to TYPE & MODEL are quoted inside ’ . . . ’ . 

Examples of the input files are shown in the figure 6.1. of 
chapter 6. 
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